Pulsating instability of a Bose-Einstein condensate in an optical lattice 
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We find numerically that in the limit of weak atom-atom interactions a Bose-Einstein condensate 
in an optical lattice may develop a pulsating dynamical instability in which the atoms nearly peri- 
odically form a peak in the occupation numbers of the lattice sites, and then return to the unstable 
initial state. Multiple peaks behaving similarly may also occur. Simple arguments show that the 
pulsating instability is a remnant of integrability, and give a handle on the relevant physical scales. 
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Wave propagation combined with nonlinearity fre- 
quently gives rise to what is known as modulational or 
dynamical instability. It has been observed experimen- 
tally pQ that a one-dimensional Bose-Einstein condensate 
(BEC) with attractive interactions between the atoms is 
unstable, and may break up into a chain of solitons |2J. 
Instabilities result from the interplay between dispersion 
and nonlinearity, and confinement of a condensate into 
an optical lattice modifies the dispersion so that the 
range of potential nonlinear phenomena is correspond- 
ingly broader [3] . The flow of a BEC in an optical lattice 
is dynamically unstable [H [5] at certain quasimomcnta 
even for repulsive interactions. The empirical view is 
that, upon the instability, the BEC in a lattice evolves 
irregularly However, one expects [7] and finds experi- 
mentally [8] that the instability may again lead to soliton 
formation. Although our focus is on a BEC in an optical 
lattice, similar phenomenology has been studied in a wide 
range of nonlinear systems in many field of science. For 
instance, there is a precise analog of the experiment [5] 
in nonlinear optics that employs a waveguide array [§] . 

In this Letter we study a BEC in an optical lattice 
numerically for weaker atom-atom interactions than is 
customary [ID] . We find that the system may exhibit a 
pulsating dynamical instability in which the atoms nearly 
periodically collect to a peak in lattice occupation num- 
bers, and subsequently disperse back to (very close to) 
the initial unstable state. Apparently related pulsations 
starting from an already compressed atom distribution 
in a lattice have been reported from numerical studies 
of the nonpolynomial Schrodinger equation [IT] , but the 
present framework affords many additional insights. By 
comparing with a two-site lattice, we see that, while the 
multisite lattice retains the instability as a vestige of the 
nonlinearity, it behaves approximately as if it were in- 
tegrable when the instability unfolds. We arc also able 
to develop simple quantitative estimates that delineate 
the conditions for the pulsating instability, and correctly 
give the temporal and spatial scales of the pulsating atom 
distribution. 

The discrete nonlinear Schrodinger equation (DNLSE) 
serves as a tight-binding model for the amplitudes of 



the BECs a*; at the sites k = 0, . . . , N — 1 of a one- 
dimensional optical lattice, 

ia k = -±5(a k +i + a k -i) + x\ a k\ 2 ctk ■ (1) 

We use periodic boundary conditions. The quantities 
|afc| 2 are proportional to the numbers of atoms at the 
sites k, and we refer to them as populations. The con- 
stants 5 and x characterize tunneling of the atoms from 
site to site and on-site atom-atom interactions. Momen- 
tarily we discuss the case with 5 > and x > 0. DNLSE 
is the equation of motion under the Hamiltonian 

H = E + afc-i) + IxKI 4 ] (2) 

k 

and the Poisson brackets {a k ,al,} = —i6k,k', so that 
the value of the Hamiltonian, the energy, is a constant 
of the motion. The normalization ^ fc \a k \ 2 of the BEC 
amplitudes, which is proportional to the total number of 
atoms, is also a constant of the motion. We normalize 
to the number of lattice sites ./V: ^ fc | ck^ | 2 = N. For 
notational simplicity N is taken to be even. 

The DNLSE admits stationary solutions of the form 
a k (p,t) = e i[pk-Kp)t]^ with n(p) = -Scosp + x- By 
virtue of the periodic boundary conditions the quasimo- 
menta must be of the form p = 2nP/N with an inte- 
ger P. In the reduced zone scheme all integers char- 
acterizing plane wave modes are folded into the inter- 
val (— N/2, N/2]. The standard ansatz for linear stabil- 
ity analysis for the plane- wave mode p is a k (t; p, q) = 
e i[pk-»(p)t}[ ue i(qk-ut) + v * e -i(qk-u*t)^ where u and v 

are (small) constants. By the periodic boundary con- 
ditions the quasimomcntum of an excitation q (displace- 
ment from the quasimomentum of the stationary solu- 
tion) must again be of the form q = 2-kQ/N with an 
integer Q ^ 0. This ansatz works, giving both the small- 
oscillation frequencies [5] 

uj{p,q) — 5 smpsvaq 
± \/<5cosp(l — cosq)[(2x + <5cosp(l — cosg)] (3) 

and the amplitudes u and v as a solution to an eigenvalue 
problem. 
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FIG. 1: Snapshots of the population \<Xk\ 2 at different times. 
The parameters are TV = 32, x/S = 0.015, and p — n. The 
instability is seeded with Gaussian noise with the amplitude 

e = io- 4 . 

The expression inside the square root in Eq. ^ is 
nonnegative for p G [~| 7r j| 7r ] but may become nega- 
tive otherwise, depending on the values of 5, \i P an 1- 
Correspondingly, the square root may become imaginary, 
whereupon there are small deviation from the steady 
state that grow exponentially in time. The flow with 
the quasimomentum p is then dynamically unstable. If 
there is an imaginary part to the eigenfrcquency lo, it is 
the same for excitation quasimomenta q and —q. More- 
over, the unstable direction, the eigenvector [u,i>] T cor- 
responding to the frequency with a negative imaginary 
part, is the same for q and —q. In this sense the modes 
q and — q are equivalent as it comes to the instability. 

In our renewed look into the fate of an unstable mode 
we carry out numerical experiments on the DNLSE. Fig- 
ure [l] shows an example result. Here we have N — 32 lat- 
tice sites and an initial flow with p — it, or P = N/2 = 16. 
In this state a k (p, t) the sign of the condensate amplitude 
simply flips from one site to the next. The interaction 
is weak, x/5 — 0.015. To speed up the instability, we 
add complex-valued Gaussian noise to the amplitude of 
each site. Here the root-mean-square noise amplitude is 
£ = 10~ 4 . The Figure shows a histogram of the popu- 
lations \dk\ 2 at certain snapshot times. The instability 
leads to a single-peaked distribution of the occupation 
numbers. Moreover, upon further time evolution the sys- 
tem returns very close to the initial state, again pulsates 
to a peak, and so on. We have periodic peaking and 
recurrences of the unstable initial state. 

Another view on the same data is given in Fig. [2j We 
plot the fraction of the initial state a k (0) in the state at 
time t, a k {t), the quantity f(t) = \J2 k a* k (0)a k (t)\ 2 /N 2 , 
as a function of time. The times of the first four deep 
minima are the times picked for the snapshots of the 
four peaks in Fig.[l] but Fig. [2] reveals five more peaking 
events. As appropriate for phenomena initiated by a dy- 
namical instability, if everything else is held unchanged, 
the time of the first peak depends logarithmically on the 




FIG. 2: Fraction of the initial state in the state of the lattice, 
/, plotted as a function of time. This figure is from the same 
data as Fig. [I] 

amplitude of the noise £. A closer inspection reveals that 
the interval between the bottoms of the dips in Fig. [2] 
varies slightly, so that the subsequent peaking events are 
not strictly periodic. 

From data sets of this kind two observations emerge. 
First, if the peaking starts close to periodic, it may con- 
tinue so seemingly indefinitely. We have examples with 
50 recurrences. Second, near-periodic recurrences occur 
also for other initial flow states than the edge of the first 
Brillouin zone with p = it. In such a case the peaks 
recur at different positions in the lattice, and move at 
approximately the usual group velocity v g = S sin p. 

The common predicate for quasiperiodic recurrences 
of a single peak is that, except for the +Q and — Q 
equivalence, one and only one of the possible modes Q 
is unstable. In the case of more than a few lattice sites, 
this condition can only be satisfied for weak interactions, 
X/S <C 1, and the one possible unstable mode is Q — 1. A 
straightforward analysis of the eigenvalues ([3| shows that 
the range of the interaction strength where the Q = 1 
mode is unstable and the Q = 2 mode is not, in the limit 
N 3> 1, is given by 



| cosp|(7r/7V) 2 < x/6 < 4 I cosp\(ir/N) 2 



(4) 



A plot similar to Fig. [T] shows that the Fourier trans- 
form of the site amplitudes a k (t) with respect to the 
site index k mainly displays the Fourier components P 
and P ± 1. This is as should be, since the excitation 
modes that go unstable have the indices Q — ±1 with 
respect to the initial plane wave and therefore make the 
Fourier components P ± 1. Depending on the parame- 
ters, the P ± 1 modes may almost completely extinguish 
the original mode P. The nonlinear dynamics also gen- 
erates Fourier component P ± 2, etc., but with smaller 
amplitudes. If all atoms were in the plane wave modes 
P± 1, we would have a modulation of the populations of 
the form cos 2 (2nk / N + ip). This neatly explains why the 
width of the base of a recurring peak in the number of 
sites k is about half of the size of the lattice N. In fact, 
let us use the lower limit of the interaction strength in 
Eq. Q as a condition for the occurrence of a peak and 
the observation that the width of the peak is about N/2, 
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then we have an estimate for the width of the peak as a 
function of the interaction strength \, 

*~7r>/<y|cosp|/x:/2. (5) 

We may then cast the condition for the occurrence of 
a single peak into the form that for a given interaction 
strength and the ensuing width of the peak £, one and 
only one peak comes up if and only if the length of the 
lattice is such that the lattice loosely accommodates one 
peak, but two peaks would have to overlap. 

Clean recurrences do not occur for small site numbers 
> 2, such as N — 4. Then a fairly strong nonlinearity 
x/S ~ 1 is required for the instability, and many Fourier 
modes are coupled in. Besides, the Fourier components 
wrap around; P + 2 and P — 2 are the same plane wave 
mode for N = 4. Apparently the nonlinearity strongly 
mixes all modes, which leads to irregular behavior. 

In search of a qualitative explanation of the quasiperi- 
odic behavior we study the system with just two sites. 

The transformation ao.i — e t ^ ± 2 A< ^'^/ ± An, with 
(n, <j>) and (An, A<p) being usual canonical-conjugate 
pairs (Poisson brackets {n, 4>} = 1 5 and so on) puts the 
two-site Hamiltonian into the form 

H 2 = x(\n 2 + An 2 ) - 8\fn 2 - 4An 2 cos A0 . (6) 

Since there is no angle <j> in the Hamiltonian, the canoni- 
cal conjugate n, total normalization, is a constant of the 
motion, and n = N = 2 according to our convention. An 
and A<f> are, respectively, half of the population difference 
and the BEC phase difference between the sites and 1. 
The potentially unstable steady state is P — 1, which in 
the present variables translates to An — 0, A(f> = it. The 
onset of the instability is at %/<5 = 1. 

The Hamiltonian ^ is also a constant of the motion, 
which makes the two-site system integrable. The mo- 
tion in phase space is along constant-energy curves in the 
(An, Acf>) plane. Accordingly, in Fig.|3]we draw constant- 
energy contours in this plane for x/S = 0.5 (a), and 
x/S = 1.5 (b). In these figures the A<j> axis runs from 
to 2tt, so that the potentially unstable steady state 
is at the center of the contour plot; the unconditionally 
stable steady state with P = 0, or An = 0, Acj) = is 
split between the middle of the lower and upper edges of 
the contour plot. For x/S = 0.5 the potentially unstable 
steady state is an elliptic fixed point, and time evolu- 
tion starting in its vicinity takes the system periodically 
around the stationary state. At x/S = I the elliptic fixed 
point bifurcates, and for x/S — 1-5 there is a homoclinic 
point plus directions along which the system either re- 
cedes from the homoclinic point exponentially in time, or 
approaches the homoclinic point according to a negative- 
exponential law. The directions of the flow are indicated 
in panel (b). Thus, starting in the vicinity of what used 
to be the potentially unstable steady state, the system 
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FIG. 3: Contour plots for the Hamiltonian in the case of two 
sites only, in the (An,A0) plane, for the interactions strengths 
(a) x/5 = 0.5 and (b) 1.5. The arrows in panel (b) indicate 
the direction of flow. 

takes off in an unstable direction and goes around one of 
the bifurcated elliptic fixed points. Depending on where 
the system starts out, it may have to loop around the 
second new elliptic fixed point, too, before returning to 
the original state. 

Multisite systems share the property of the two-site 
system that, starting from noise in the neighborhood of 
the unstable steady state, they may execute a large ex- 
cursion out of the vicinity of the steady state but then 
return. We surmise that the mechanism of the recur- 
rences in multisite and two-site systems is basically the 
same. As the interaction strength increases, the would- 
be unstable fixed point bifurcates and generates a ho- 
moclinic point along with two elliptic fixed points; ac- 
tually, two overlapping homoclinic points for the modes 
Q = ±1. Starting from noise, the system is propelled 
along a newly created unstable direction, which is the 
same for both modes Q — ±1. A trajectory starting 
from the homoclinic point in an unstable direction re- 
turns to the homoclinic point along a stable direction. 
The multisite system appears to stay in the vicinity of 
such a homoclinic orbit. 

Depending on the motion in the neighborhood of the 
homoclinic orbit due to the initial noise, the multisite 
system does not return to exactly where it started from. 
This may account for some of the variations in the period. 
Nonetheless, for small initial noise, the system should 
spend most of its time either receding from the homo- 
clinic point, or approaching it. Accordingly, we have ver- 
ified by numerical experiments that the time scale of the 
instability r = l/|Qxt;| is the dominant time scale for the 
recurrences. 

Although the two-site model offers a neat overall ex- 
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FIG. 4: (color online). Density plot for the populations \ctk\ 2 
as a function of time for the parameters N = 128, x/<5 = 
0.015, and £ = 10 -4 ; the darker the area, the smaller the 
value. 

planation of the pulsating instability, theory of nonlinear 
systems should offer a more fine-graded description. We 
could have here a g-breather [12] , or an oscillatory insta- 
bility P2J H3] that takes the lattice between the unsta- 
ble plane wave and a breather. The difference between 
the time for the first peak and the intervals between the 
subsequent peaks could be because the system initially 
radiates the energy that does not participate in the os- 
cillations into the background plane wave modes. The 
details remain to be investigated. 

For a lattice with a large number of sites it may be diffi- 
cult to realize the one-peak condition Q experimentally, 
but there is some leeway: the noise may seed several ap- 
proximately independent pulsating peaks in the lattice. 
Figure [4]shows an example. Here we have N = 128 lattice 
sites, the interaction strength is x/$ = 0.015, and initial 
random noise on the p = it mode is specified by £ = 10 -4 . 
The populations \cek\ 2 are presented on a density plot as 
a function of time, the darkest color indicating |ckfc| 2 ~ 
and white the largest populations. The right edge of the 
plot wraps around to the left edge by virtue of the peri- 
odic boundary conditions. The peaks represented by the 
white spots are not independent, they move around, join, 
and split, but most of the time there clearly are four pul- 
sating peaks. For the parameters of Fig. [4] the width of 
the peaks from Eq. ^ is £ ~ 13, so that four peaks and 
a peak's width of average spacing between them would 
require about 100 lattice sites, five peaks 130. The peak- 
width argument agrees with the numerical observations. 

Flow states with p ~ ir might be prepared, e.g., by ac- 
celerating the lattice [8], and there is also an alternative 
way. Namely, if ak(t) are the solution to the equation 
of motion (nj for the parameters (5,x), then (— \) k a* k {t) 



are the solution for the parameters (S, — x). The ground 
state p — for a lattice with repulsive interactions be- 
haves like the state p = ir after the sign of the interac- 
tions is suddenly flipped. This type of an experiment has 
been carried out using a Feshbach resonance [T], albeit 
in a nonlattice gas and with strong interactions, giving 
rise to a chain of stable solitons. However, the pulsat- 
ing instability occurs at weak interactions, and the cor- 
responding time scales may present severe technical and 
fundamental challenges. For instance, quantum fluctu- 
ations [TS] may have to be considered. Given that the 
pulsations are a result of near-integrability, we speculate 
that similar phenomena take place in lattice with hard- 
wall boundary conditions and even in the presence of soft 
trapping in the direction of the lattice. 

There is more to the instability of an unstable flow in 
an optical lattice than devolution into irregular behavior 
or soliton formation. In the case when the nonlinearity is 
weak, the instability may lead to quasiperiodic pulsation 
where the atoms first gather and then return to the ini- 
tial unstable state. With an increasing number of lattice 
sites, one may also see multiple peaks in atom popula- 
tions oscillating in a similar way. A picture based on the 
properties of the phase space of a two-site lattice qualita- 
tively explains the quasiperiodic revivals as a remnant of 
the integrability of the two-site system. Exceedingly sim- 
ple arguments give a semi-quantitative description useful 
for the assessment of the feasibility, and for the design, 
of real experiments. 
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